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Abstract. We consider the approximation properties of finite element spaces on quadri- 
lateral meshes. The finite element spaces are constructed starting with a given finite di- 
mensional space of functions on a square reference element, which is then transformed to a 
space of functions on each convex quadrilateral clement via a bilinear isomorphism of the 
square onto the element. It is known that for afhne isomorphisms, a necessary and suffi- 
cient condition for approximation of order r + 1 in L 2 and order r in H 1 is that the given 
space of functions on the reference element contain all polynomial functions of total degree 
at most r. In the case of bilinear isomorphisms, it is known that the same estimates hold 
if the function space contains all polynomial functions of separate degree r. We show, by 
means of a counterexample, that this latter condition is also necessary. As applications we 
demonstrate degradation of the convergence order on quadrilateral meshes as compared to 
rectangular meshes for serendipity finite elements and for various mixed and nonconforming 
finite elements. 



1. Introduction 

Finite element spaces are often constructed starting with a finite dimensional space V of 
shape functions given on a reference element K and a class S of isomorphic mappings of the 
reference element. If F G S we obtain a space of functions Vp(K) on the image element 
K = F(K) as the compositions of functions in V with F" 1 . Then, given a partition 7 of 
a domain Q into images of K under mappings in S, we obtain a finite element space as a 
subspaceQ of the space V 7 of all functions on Q which restrict to an element of Vf(K) on 
each K G 7. 

For example, if the reference element K is the unit triangle, and the reference space V 
is the space V r (K) of polynomials of degree at most r on K, and the mapping class S is 
the space Aff(i\~) of affine isomorphisms of K into IR 2 , then V 7 is the familiar space of 
all piecewise polynomials of degree at most r on an arbitrary triangular mesh 7. When 
S = Aff(i\~), as in this case, we speak of affine finite elements. 

If the reference element K is the unit square, then it is often useful to take S equal to a 
larger space than Aff(i\~), namely the space Bil(i\~) of all bilinear isomorphisms of K into IR 2 . 
Indeed, if we allow only affine images of the unit square, then we obtain only parallelograms, 
and we are quite limited as to the domains that we can mesh (e.g., it is not possible to mesh 
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1 The subspace is typically determined by some interelement continuity conditions. The imposition of such 
conditions through the association of local degrees of freedom is an important part of the construction of 
finite element spaces, but, not being directly relevant to the present work, will not be discussed. 
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a triangle with parallelograms). On the other hand, with bilinear images of the square we 
obtain arbitrary convex quadrilaterals, which can be used to mesh arbitrary polygons. 

The above framework is also well suited to studying the approximation properties of finite 
element spaces. See, e.g., [0 and A fundamental result holds in the case of affine finite 
elements: S = Aff(A'). Under the assumption that the reference space V D V r {K), the 
following result is well known: if Ti, T2, ... is any shape-regular sequence of triangulations 
of a domain Q and u is any smooth function on Q, then the L 2 error in the best approximation 
of u by functions in V 7 " is 0(h r+1 ) and the piecewise H l error is 0(h r ), where h = h{7 n ) is 
the maximum element diameter. It is also true, even if less well-known, that the condition 
that V D V r (K) is necessary if these estimates are to hold. 

The above result does not restrict the choice of reference element K, so it applies to 
rectangular and parallelogram meshes by taking K to be the unit square. But it does not 
apply to general quadrilateral meshes, since to obtain them we must choose S = BW(K), 
and the result only applies to affine finite elements. In this case there is a standard result 
analogous to the positive result in the previous paragraph, [0, Jj]], |4], Section I. A. 2]. Namely, 
if V D Q r (K), then for any shape-regular sequence of quadrilateral partitions of a domain Q 
and any smooth function u on Q, we again obtain that the error in the best approximation 
of u by functions in V 7 " is 0(h r+1 ) in L 2 and 0{h r ) in (piecewise) H 1 . It turns out, as 
we shall show in this paper, that the hypothesis that V D Q r {K) is strictly necessary for 
these estimates to hold. In particular, if V D V r {K) but V ~j) Q r (K), then the rate of 
approximation achieved on general shape-regular quadrilateral meshes will be strictly lower 
than is obtained using meshes of rectangles or parallelograms. 

More precisely, we shall exhibit in Section 3 a domain Q and a sequence, Ti, T2, ... of 
quadrilateral meshes of it, and prove that whenever V(K) ~jt) Q r (K), then there is a function 
u on Q such that 

inf ||u - v\\ L 2 (n) ^ o{h r ), 

(and so, a fortiori, is 7^ 0(h r+1 )). A similar result holds for H 1 approximation. The 
counterexample is far from pathological. Indeed, the domain Q is as simple as possible, 
namely a square; the mesh sequence T n is simple and as shape-regular as possible in that all 
elements at all mesh levels are similar to a single fixed trapezoid; and the function u is as 
smooth as possible, namely a polynomial of degree r. 

The use of a reference space which contains V r (K) but not Q r (K) is not unusual, but 
the degradation of convergence order that this implies on general quadrilateral meshes in 
comparison to rectangular (or parallelogram) meshes is not widely appreciated. It has been 
observed in special cases, often as a result of numerical experiments, cf. JF], Section 8.7]. 

We finish this introduction by considering some examples. Henceforth we shall always 
use K to denote the unit square. First, consider finite elements with the simple polynomial 
spaces as shape functions: V = V r (K). These of course yield 0(h r+1 ) approximation in L 2 
for rectangular meshes. However, since V T [K) D Q[ r /2](K) but V r (K) ^ Q|_ r /2j+i(-^), on 
general quadrilateral meshes they only afford 0(/i' r / 2 ' +1 ) approximation. 

A similar situation holds for serendipity finite element spaces, which have been popular 
in engineering computation for thirty years. These spaces are constructed using as reference 
shape functions the space S r (K) which is the span of V r {K) together with the two monomials 
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x r y and yx r . (The purpose of the additional two functions is to allow local degrees of freedom 
which can be used to ensure interelement continuity.) For r = 1, S\(K) = Q\(K), but for 
r > 1 the situation is similar to that for V r (K), namely S r (K) Q^ r / 2 \(K) but S r (K) ^ 
Q\r/2\+i(K). So, again, the asymptotic accuracy achieved for general quadrilateral meshes 
is only about half that achieved for rectangular meshes: 0(/i' r / 2 l +1 ) in L 2 and 0(/iL r / 2 J) in 
H l . In Section 4 we illustrate this with a numerical example. 

While the serendipity elements are commonly used for solving second order differential 
equations, the pure polynomial spaces V r can only be used on quadrilaterals when interele- 
ment continuity is not required. This is the case in several mixed methods. For example, a 
popular element choice to solve the stationary Stokes equations is bilinearly mapped piece- 
wise continuous <2 2 elements for the two components of velocity, and discontinuous piecewise 
linear elements for the pressure. Typically the pressure space is taken to be functions which 
belong to Vi(K) on each element K. This is known to be a stable mixed method and gives 
second order convergence in H l for the velocity and L 2 for the pressure. If one were to define 
the pressure space instead by using the construction discussed above, namely by composing 
linear functions on reference square with bilinear mappings, then the approximation prop- 
erties of mapped V\ discussed above would imply that method could be at most first order 
accurate, at least for the pressures. Hence, although the use of mapped V\ as an alternative 
to unmapped V\ pressure elements is sometimes proposed 0, it is probably not advisable. 

Another place where mapped V r spaces arise is for approximating the scalar variable in 
mixed finite element methods for second order elliptic equations. Although the scalar variable 
is discontinuous, in order to prove stability it is generally necessary to define the space for 
approximating it by composition with the mapping to the reference element (while the space 
for the vector variable is defined by a contravariant mapping associated with the mapping 
to the reference element). In the case of the Raviart-Thomas rectangular elements, the 
scalar space on the reference square is Q r (K), which maintains full 0(h r+1 ) approximation 
properties under bilinear mappings. By contrast, the scalar space used with the Brezzi- 
Douglas-Marini and the Brezzi-Douglas-Fortin-Marini spaces is V r (K). This necessarily 
results in a loss of approximation order when mapped to quadrilaterals by bilinear mappings. 

Another type of element which shares this difficulty is the simplest nonconforming quadri- 
lateral element, which generalizes to quadrilaterals the well-known piecewise linear non- 
conforming element on triangles, with degrees of freedom at the midpoints of edges. On the 
square, a bilinear function is not well-defined by giving its value at the midpoint of edges (or 
its average on edges), because these quantities do not comprise a unisolvent set of degrees of 
freedom (the function (x — 1/2) (y — 1/2) vanishes at the four midpoints of the edges of the 
unit square). Hence, various definitions of nonconforming elements on rectangles replace the 
basis function xy by some other function such as x 2 —y 2 . Consequently, the reference space 
contains Vi(K), but does not contain Q 1 (K), and so there is a degradation of convergence 
on quadrilateral meshes. This is discussed and analyzed in the context of the Stokes problem 
in @. 

As a final application, we remark that many of the finite element methods proposed for 
the Reissner-Mindlin plate problem are based on mixed methods for the Stokes equations 
and/or for second order elliptic problems. As a result, many of them suffer from the same 
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sort of degradation of convergence on quadrilateral meshes. An analysis of a variety of these 
elements will appear in forthcoming work by the present authors. 

In Section 3, we prove our main result, the necessity of the condition that the reference 
space contain Q r (K) in order to obtain 0(h r+1 ) approximation on quadrilateral meshes. The 
proof relies on an analogous result for affine approximation on rectangular meshes, where 
the space V r (K) enters rather than Q r (K). While this is a special case of known results, 
for the convenience of the reader we include an elementary proof in Section 2. In the final 
section we illustrate the results with numerical computations. 

2. Approximation theory of rectangular elements 

In this section we prove some results concerning approximation by rectangular elements 
which will be needed in the next section where the main results are proved. The results in 
this section are essentially known, and many are true in far greater generality than stated 
here. 

If K is any square with edges parallel to the axes, then K = F K (K) where F K (x) : = 
xk + hxx with xk G and hx > the side length. For any function u G L 2 (K), we define 
u K = u o Fk G L 2 [K\ i.e., uk(x) = u(xk + hxx). Given a subspace S of L 2 (K) we define 
the associated subspace on an arbitrary square K by 

S(K) = {u: K ^R\u K e S}. 

Finally, let Q denote the unit cube (Q and K both denote the unit cube, but we use the 
notation Q when we think of it as a fixed domain, while we use K when we think of it as 
a reference element). For n — 1, 2, . . . , let be the uniform mesh of Q into n d subcubes 
when h = 1/n, and define 

S h = { u : Q -> R | u\ K e S(K) for all K e7 h }. 

In this definition, when we write u\k £ S{K) we mean only that u\k agrees with a function 
in Sk almost everywhere, and so do not impose any interelement continuity. 

The following theorem gives a set of equivalent conditions for optimal order approximation 
of a smooth function u by elements of Sh- 

Theorem 1. Let S be a finite dimensional subspace of L 2 (K), r a non-negative integer. 
The following conditions are equivalent: 

1. There is a constant C such that inf ||u — i>||i,2(m < Ch r+1 \u\ r+ i for all u e H r+1 (Q). 

2. inf ||w — v \\l 2 (q.) = °{h r ) for all u G V r (£l). 

ves h 

3. V r (K) C S. 

Proof. The first condition implies that mf ve s h \\u — v\\ L 2^ = for u G V r (Q), and so implies 
the second condition. The fact that the third condition implies the first is a well-known 
consequence of the Bramble-Hilbert lemma. So we need only show that the second condition 
implies the third. 

The proof is by induction on r. First consider the case r = 0. We have 
(1) inf \\u - u||£ 2( n) = Yl [n LJ u ~ Vk Wl 2 (k) = h 2 inf JI^ ~ Hl^tfy 
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where we have made the change of variable w = vr in the last step. 
In particular, for u = 1 on Q, uk = 1 on K for all K, so the quantity 



is independent of K. Thus 



c := inf \\u K - w\\ 

wes { ' 



jnf \\u-v\\ 2 2{Q) = h 2 J2 c 

Kei h 



The hypothesis that this quantity is o(l) implies that c = 0, i.e., that the constant function 
belongs to S. 

Now we consider the case r > 0. We again apply ([!]), this time for u an arbitrary 
homogeneous polynomial of degree r. Then 

(2) Uk{x) = u(xk + hx) = u{hx) + p(x) = h r u(x) +p(x), 

where p G V r -\{K\ Substituting in (ffl), and invoking the inductive hypothesis that S D 
Vr-i(K), we get that 

i?J \\u - v\\ 2 L2{n) = h 2+2r inf J^ - w \\% ( k) = h2r inf J U - w W 2 mky 

KEJh 

Again the last infimum is independent of K so we immediately deduce that u belongs to S. 
Thus S contains all homogeneous polynomials of degree r and all polynomials of degree less 
than r (by induction), so it indeed contains all polynomials of degree at most r. □ 

A similar theorem holds for gradient approximation. Since the finite elements are not 
necessarily continuous we write for the gradient operator applied piecewise on each 
element. 

Theorem 2. Let S be a finite dimensional subspace of L 2 (K), r a non-negative integer. 
The following conditions are equivalent: 

1. There is a constant C such that inf HV/^w — f)IU 2 (n) — Ch r \u\ r+ \ for all u G H r+1 (Q). 

v<=S h 

2. inf ||V/t(u — w)||l2(q) = o(h r ~ l ) for all u G V r (VL). 

ves h 

3. V r (K) CT (K) + S. 

Proof. Again, we need only prove that the second condition implies the third. In analogy to 
([[]), we have 

mf J2\\V(u-v)\\ 2 LHK) = £^ inf \Mu-v K )\\ 2 LHK) 
1 Ker h Kei h K 

= Yl ^tW^^K -W)\\ 2 L2 £ 

where we have made the change of variable w = vk in the last step. 

The proof proceeds by induction on r, the case r = being trivial. For r > 0, apply 
(0) with u an arbitrary homogeneous polynomial of degree r. Substituting (FJ) in ([3]), and 
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invoking the inductive hypothesis that Vq(K) + S 3 V r -i(K), we get that 

inf \\V h (u-v)\\l, m = h 2r £ inf ||V(t* - w)\\ 2 = h 2r ' 2 inf ||V(u - w)\\ 2 

KeJh 

Since we assume that this quantity is o(h 2r ~ 2 ), the last infimum must be 0, so u differs 
from an element S by a constant. Thus Vq{K) + S contains all homogeneous polynomials of 
degree r and all polynomials of degree less than r (by induction), so it indeed contains all 
polynomials of degree at most r. □ 

Remarks. 1. If S contains Vo(K), which is usually the case, then the third condition of 
Theorem reduces to that of Theorem [I]. 

2. A similar result holds for higher derivatives (replace by V™ in the first two conditions, 
and V (K) by V m -i(K) in the third). 



3. Approximation theory of quadrilateral elements 

In this, the main section of the paper, we consider the approximation properties of finite 
element spaces defined with respect to quadrilateral meshes using bilinear mappings starting 
from a given finite dimensional space of polynomials V on the unit square K = [0, 1] x [0, 1]. 
For simplicity we assume that V 3 Vo(K). For example V might be the space V r (K) of 
polynomials of degree at most r, or the space Q r (K) of polynomials of degree at most r in 
each variable separately, or the serendipity space S r (K) spanned by V r (K) together with the 
monomials x\x2 and x\x r 2 . Let F be a bilinear isomorphism of K onto a convex quadrilateral 
K = F(K). Then for u G L 2 (K) we define iip G L 2 (K) by up{x) = u(Fx), and set 

V F (K) = {u: K ^R\u F eV}. 

(Note that the definition of this space depends on the particular choice of the bilinear iso- 
morphism F of K onto K, but whenever the space V is invariant under the symmetries of 
the square, which is usually the case in practice, this will not be so.) We also note that the 
functions in Vp(K) need not be polynomials if F is not affine, i.e., if K is not a parallelogram. 

Given a quadrilateral mesh T of some domain, Q, we can then construct the space of 
functions V 7 consisting of functions on the domain which when restricted to a quadrilateral 
K G T belong to Vf k {K) where Fk is a bilinear isomorphism of K onto K. (Again, if V is 
not invariant under the symmetries of the square, the space V 7 will depend on the specific 
choice of the maps F K .) 

It follows from the results of the previous section that if we consider the sequence of 
meshes of the unit square into congruent subsquares of side length h = 1/n, then each of 
the approximation estimates 

(4) inf ||u - v\\ L 2 {n) < Ch r+1 \u\ r+1 for all u G H r+ \tt), 

vev 7 h 

(5) inf \\V h (u - v)\\ L 2 {n) < Ch r \u\ r+1 for all u G H r+l (tt), 

vev' T h 

holds if and only V r (K) C V. It is not hard to extend these estimates to shape- regular 
sequences of parallelogram meshes as well. However, in this section we show that for these 
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estimates to hold for more general quadrilateral mesh sequences, a stronger condition on V 
is required, namely that V D Q r (K). 

The positive result, that when V D Q r (K), then the estimates @ and fl5|) hold for any 
shape regular sequence of quadrilateral meshes 1^, is known. See, e.g., [0, Q, or @, Section 
I. A. 2]. We wish to show the necessity of the condition V D Q r (K). 

As a first step we show that the condition Vp(K) D V r (K) is necessary and sufficient to 
have that V D Q r (K) whenever F is a bilinear isomorphism of K onto a convex quadrilateral. 
This is proven in the following two theorems. 

Theorem 3. Suppose that V D Q r (K). Let F be any bilinear isomorphism of K onto a 
convex quadrilateral. Then Vf{K) D V r {K). 

Proof. The components of F(x, y) are linear functions of x and y, so if p is a polynomial of 
degree at most r, then p(F(x, y)) is of degree at most r in x and y, i.e., p o F £ Q r (K) C V. 
Therefore p £ Vp(K). □ 

The reverse implication holds even under the weaker assumption that Vf{K) contains 
V T {K) just for the two specific bilinear isomorphism 

F(x, y) = (x, y(x + 1)), F(x, y) = (y, x(y + 1)), 

both of which map K isomorphically onto the quadrilateral K' with vertices (0,0), (1,0), 
(0, 1), and (1,2). This fact is established below. 

Theorem 4. Let V be a vectorspace of functions on K. Suppose that Q r (K) <£. V. Then 
either Vp(K') £ V r {K') or Vp(K') £ P r (K'). 

Remark. If the space V is invariant under the symmetries of the square, then Vp(K') = 
Vp(K') so neither contains V r (K f ). 

Proof. Assume to the contrary that Vp(K') D V r {K') and Vp(K') D V r (K'). We prove 
that V D Q r (K) by induction on r. The case r = being true by assumption, we consider 
r > and show that the monomials x r y s and x s y r belong to V for s = 0, 1, . . . , r. From the 
identity 

(6) x r y s = x r - s [y(x + l)] s - ^ ("V^ = vY^m, y) s - ^ (j)^ 5 , 

t=i ^ ' t=i \ / 

we see that for < s < r, the monomial x r y s is the sum of a polynomial which clearly 
belongs to V (since Fi(x,y) r ~ s F2(x,y) s = x r ~ s y s £ V r (K') C Vp(K')) and a polynomial in 
Q r -i(K), which belongs to V by induction. Thus each of the monomials x r y s with < s < r 
belongs to V, and, using F, we similarly see that all the monomials x s y r , < s < r belong 
to V. Finally, from @ with s = r, we see that x r y r is a linear combination of an element of 
V and monomials x s y r with s < r, so it too belongs to V. □ 

We now combine this result with the those of the previous section to show the necessity 
of the condition V 3 Q r (K) for optimal order approximation. Let V be some fixed finite 
dimensional subspace of L 2 (K) which does not include Q r (K). Consider the specific division 
of the unit square K into four quadrilaterals shown on the left in Figure |l|. For definiteness we 
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place the vertices of the quadrilaterals at (0, 1/3), (1/2, 2/3) and (1, 1/3) and the midpoints 
of the horizontal edges and the corners of K. 



Figure 1. a. A partition of the square into four trapezoids, b. A mesh 
composed of translated dilates of this partition. 

The four quadrilaterals are mutually congruent and affinely related to the specific quadri- 
lateral K' defined above. Therefore, by Theorem ^, we can define for each of the four 
quadrilaterals K" shown in Figure |I] an isomorphism F" from the unit square so that 
Vf"{K") ^ V r {K"). If we let 5" be the subspace of L 2 {K) consisting of functions which 
restrict to elements of Vf"{K") on each of the four quadrilaterals K", then certainly S does 
not contain V r {K), since even its restriction to any one of the quadrilaterals K" does not 
contain P r {K"). 

Next, for n = 1,2,... consider the mesh 7' h of the unit square Q shown in Figure [I]b, 
obtained by first dividing it into a uniform n x n mesh of subsquares, n = 1/h, and then 
dividing each subsquare as in Figure 0a. Then the space of functions u on Q whose restrictions 
on each subsquare K £ 7h satisfy ur(£) = u{xk + hx) with uk £ S is precisely the same as 
the space V(7' h ) constructed from the initial space V and the mesh 7' h . In view of Theorems [l] 
and ^ and the fact that S ^ V r (K), the estimates (f|) and (IP do not hold. In fact, neither 
of the estimates 

inf \\u-v\\ L 2 {n) = o(h r ), 

veV(T h ) 

nor 

M W(u-v)\\v {p) = o(hr 1 ), 
vev(i h ) 

holds, even for u £ V r (Q). 

While the condition V D Q r (K) is necessary for 0(h r+1 ) on general quadrilateral meshes, 
the conditions V ~D V r (K) suffices for meshes of parallelograms. Naturally, the same is 
true for meshes whose elements are sufficiently close to parallelograms. We conclude this 
section with a precise statement of this result and a sketch of the proof. If V D V r (K) and 
K = F(K) with F £ Bil(i^), then by standard arguments, as in M, we get 

\\v - rrv\\ L2{K) < CppW^l^v o F\ Hr+1{ky 

where Jp is the Jacobian determinant of F. Now, using the formula for the derivative of a 
composition (as in, e.g., p. 222]), and the fact that F is quadratic, and so its third and 
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higher derivatives vanish, we get that 

L(r+i)/2j 

\ VoF \ H r+^(K)< C \\ J F-A\ 1 L^ K )\\ V \\H r + 1 (K) \ F \wl~(l)\ F \wl(K)- 

Now, 

ll^llx 00 ^) — ^^Ki II Jf^ 1 II L°°(A") — Ctlp?, \F\w^{K) — CflK, 

where hx is the diameter of K and C depends only on the shape-regularity of K. We thus 
get 

\\v - ttv\\ L 2 (k) < C\\v\\ H r + i (K) hT K l ~ 2i \ F ^wUKY 

i 

It follows that if \F\ W2 r%\ = O^ti^), we get the desired estimate 

||?; - nv\\ L 2 {K) < Ch^WvWnr+HK)- 

Following [g], we measure the deviation of a quadrilateral from a parallelogram, by the 
quantity ok '■= max(|7r — 9i\, \n — 9 2 \), where 6\ is the angle between the outward normals 
of two opposite sides of K and 82 is the angle between the outward normals of the other two 
sides. Thus < ok < 7T, with ok = if and only if K is a parallelogram. As pointed out in 
II > I-^Ivk 2 (K) — ChxitiK + ok)- This motivates the definition that a family of quadrilateral 
meshes is asymptotically parallelogram if o K = 0(h K ), i.e., if ax/hx is uniformly bounded 
for all the elements in all the meshes. From the foregoing considerations, if the reference 
space contains V r {K) we obtain 0{h r+1 ) convergence for asymptotically parallelogram, shape 
regular meshes. 

As a final note, we remark that any polygon can be meshed by an asymptotically paral- 
lelogram, shape regular family of meshes with mesh size tending to zero. Indeed, if we begin 
with any mesh of convex quadrilaterals, and refine it by dividing each quadrilateral in four 
by connecting the midpoints of the opposite edges, and continue in this fashion, as in the 
last row of Figure [| the resulting mesh is asymptotically parallelogram and shape regular. 

4. Numerical results 

In this section we report on results from a numerical study of the behavior of piecewise 
continuous mapped biquadratic and serendipity finite elements on quadrilateral meshes (i.e., 
the finite element spaces are constructed starting from the spaces Q^{K) and S 2 (K) on the 
reference square, and then imposing continuity). We present the results of two test problems. 
In both we solve the Dirichlet problem for Poisson's equation 

(7) —Am = / in Q, u = g on dfl, 

where the domain Q is the unit square. In the first problem, / and g are taken so that the 
exact solution is the quartic polynomial 

u(x, y) = x 3 + 5y 2 — 10y 3 + y 4 . 

Table [l| shows results for both types of elements using meshes from each of the first two mesh 
sequences shown in Figure 0. The first sequence of meshes consists of uniform square subdi- 
visions of the domain into n x n subsquares, n — 2, 4, 8, Meshes in the second sequence 
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are partitions of the domain into n x n congruent trapezoids, all similar to the trapezoid 
with vertices (0,0), (1/2,0), (1/2,2/3), and (0,1/3). In Table we report the errors in L 2 
for the finite element solution and its gradient both in absolute terms and as a percentage of 
the L 2 norm of the exact solution and its gradient, and we also report the apparent rate of 
convergence based on consecutive meshes in a sequence. For this test problem, the rates of 
convergence are very clear: for either mesh sequence the mapped biquadratic elements con- 
verge with the expected order 3 for the solution and 2 for its gradient. The same is true for 
the serendipity elements on the square meshes, but, as predicted by the theory given above, 
for the trapezoidal mesh sequence the order of convergence for the serendipity elements is 
reduced by 1 both for the solution and its gradient. 




FIGURE 2. Three sequences of meshes of the unit square: square, trapezoidal, and 
asymptotically parallelogram. Each is shown for n — 2 ,4, 8, and 16. 

As a second test example we again solved the Dirichlet problem (0), but this time choosing 
the data so that the solution is the sharply peaked function 

u(x, y) = exp(-100[(x - 1/4) 2 + (y - 1/3) 2 ]) . 

As seen in Table f|, in this case the loss of convergence order for the serendipity elements on 
the trapezoidal mesh is not nearly as clear. Some loss is evident, but apparently very fine 
meshes (and very high precision computation) would be required to see the final asymptotic 
orders. 

Finally we return to the first test problem, and consider the behavior of the serendipity 
elements on the third mesh sequence shown in Figure 0. This mesh sequence begins with 
the same mesh of four quadrilaterals as in previous case, and continues with systematic 
refinement as described at the end of the last section, and so is asymptotically parallelogram. 
Therefore, as explained there, the rate of convergence for serendipity elements is the same 
as for affine meshes. This is clearly illustrated in Table []. 

While the asymptotic rates predicted by the theory are confirmed in these examples, it is 
worth noting that in absolute terms the effect of the degraded convergence rate is not very 
pronounced. For the first example, on a moderately fine mesh of 16 x 16 trapezoids, the 
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TABLE 1. Errors and rates of convergence for the test problem with polyno- 
mial solution. 



Mapped biquadratic elements 
square meshes trapezoidal meshes 







||u- 


U h \\ L 2 




l|V(« 


- U h )\\ L 2 






||u- 






||V(« 


- U h )\\ L 2 




n 


err. 


% 


rate 


err. 


% 


rate 


err. 


% 


rate 


err. 


% 


rate 


2 


3.5e- 


02 


2.877 




4.5e-01 


37.253 




4.8e- 


02 


3.951 




5.9e-01 


48.576 




4 


4.4e- 


03 


0.360 


3.0 


l.le-01 


9.333 


2.0 


5.8e- 


03 


0.475 


3.1 


1.5e-01 


12.082 


2.0 


8 


5.5e- 


04 


0.045 


3.0 


2.8e-02 


2.329 


2.0 


7.1e- 


04 


0.058 


3.0 


3.7e-02 


3.017 


2.0 


16 


6.9e- 


05 


0.006 


3.0 


7.1e-03 


0.583 


2.0 


8.7e- 


05 


0.007 


3.0 


9.2e-03 


0.753 


2.0 


32 


8.6e- 


06 


0.001 


3.0 


1.8c-03 


0.146 


2.0 


1.1c- 


05 


0.001 


3.0 


2.3e-03 


0.188 


2.0 


64 


l.le- 


06 


0.000 


3.0 


4.4e-04 


0.036 


2.0 


1.3e- 


06 


0.000 


3.0 


5.7e-04 


0.047 


2.0 



Serendipity elements 
square meshes trapezoidal meshes 







ll«- 


uh\\l* 




l|V(« 


- U h )\\ L 2 






ll«- 


UhWh* 




l|V(« 


- U h )\\ L 2 




n 


err. 


% 


rate 


err. 


% 


rate 


err. 


% 


rate 


err. 


% 


rate 


2 


3.5e- 


02 


2.877 




4.5e-01 


37.252 




5.0e- 


02 


4.066 




6.2e-01 


51.214 




4 


4.4e- 


03 


0.360 


3.0 


l.le-01 


9.333 


2.0 


6.7e- 


03 


0.548 


2.9 


1.8e-01 


14.718 


1.8 


8 


5.5e- 


04 


0.045 


3.0 


2.8e-02 


2.329 


2.0 


9.7e- 


04 


0.080 


2.8 


5.9e-02 


4.836 


1.6 


16 


6.9e- 


05 


0.006 


3.0 


7.1e-03 


0.583 


2.0 


1.6e- 


04 


0.013 


2.6 


2.3e-02 


1.890 


1.4 


32 


8.6e- 


06 


0.001 


3.0 


1.8c-03 


0.146 


2.0 


3.3e- 


05 


0.003 


2.3 


1.0e-02 


0.842 


1.2 


64 


l.le- 


06 


0.000 


3.0 


4.4e-04 


0.036 


2.0 


7.4e- 


06 


0.001 


2.1 


4.9e-03 


0.401 


1.1 



solution error with serendipity elements exceeds that of mapped biquadratic elements by a 
factor of about 2, and the gradient error by a factor of 2.5. Even on the finest mesh shown, 
with 64 x 64 elements, the factors are only about 5.5 and 8.5, respectively. Of course, if 
we were to compute on finer and finer meshes with sufficiently high precision, these factors 
would tend to infinity. Indeed, on any quadrilateral mesh which contains a non-parallelogram 
element, the analogous factors can be made as large as desired by choosing a problem in which 
the exact solution is sufficiently close to — or even equal to — a quadratic function, which the 
mapped biquadratic elements capture exactly, while the serendipity elements do not (such a 
quadratic function always exists). However, it is not unusual that the serendipity elements 
perform almost as well as the mapped biquadratic elements for reasonable, and even for 
quite small, levels of error. This, together with their optimal convergence on asymptotically 
parallelogram meshes, provides an explanation of why the lower rates of convergence have 
not been widely noted. 
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TABLE 2. Errors and rates of convergence for the test problem with expo- 
nential solution. 

Mapped biquadratic elements 



square meshes trapezoidal meshes 







|| U 


- U h \\ L 2 




l|V(« 


- U h )\\ L 2 






|| U 


- U h \\ L 2 




l|V(« 


- U h )\\ L 2 




n 


err. 


% 


rate 


err. 


% 


rate 


err. 


% 


rate 


err. 


% 


rate 


2 


2.8e- 


01 


224.000 




3.0e+00 


169.630 




2.6e- 


01 


204.800 




2.8c+00 


159.208 




4 


1.2e- 


01 


93.600 


1.3 


1.5c+00 


87.322 


1.0 


2.1e- 


01 


169.600 


0.3 


1.8e+00 


99.305 


0.7 


8 


1.7e- 


02 


13.520 


2.8 


4.6e-01 


25.809 


1.8 


2.3e- 


02 


18.160 


3.2 


5.9e-01 


33.185 


1.6 


16 


1.1c- 


03 


0.920 


3.9 


1.0e-01 


5.860 


2.1 


1.3e- 


03 


1.048 


4.1 


1.2e-01 


6.819 


2.3 


32 


1.3e- 


04 


0.101 


3.2 


2.5e-02 


1.424 


2.0 


1.5e- 


04 


0.124 


3.1 


3.2e-02 


1.794 


1.9 


64 


1.5e- 


05 


0.012 


3.1 


6.3e-03 


0.354 


2.0 


1.9e- 


05 


0.015 


3.0 


7.9e-03 


0.448 


2.0 


128 


1.9e- 


06 


0.002 


3.0 


1.6e-03 


0.088 


2.0 


2.4e- 


06 


0.002 


3.0 


2.0c-03 


0.112 


2.0 



Serendipity elements 
square meshes trapezoidal meshes 







\\u 


- U h \\ L 2 




1 


V(t/ 


- U h )\\ L 2 






\\u 


- U h \\ L 2 




1 


V(u 


- U h )\\ L 2 




n 


err 




% 


rate 


err 




% 


rate 


err 




% 


rate 


err 




% 


rate 


2 


2.0e- 


01 


159.200 




2.4c4 


-00 


133.372 




2.1e- 


01 


169.600 




2.3c4 


-00 


130.340 




4 


1.2e- 


01 


92.000 


0.8 


1.4eH 


-00 


80.531 


0.7 


2.1e- 


01 


168.000 


0.0 


1.7e4 


-00 


93.819 


0.5 


8 


1.7e- 


02 


13.520 


2.8 


4.6e- 


-01 


26.293 


1.6 


2.4e- 


02 


18.880 


3.2 


6.1e- 


-01 


34.564 


1.4 


16 


l.le- 


03 


0.920 


3.9 


l.le- 


-01 


5.948 


2.1 


1.5e- 


03 


1.208 


4.0 


1.4e- 


-01 


7.737 


2.2 


32 


1.3e- 


04 


0.101 


3.2 


2.5e- 


-02 


1.432 


2.1 


2.0e- 


04 


0.162 


2.9 


3.8e- 


-02 


2.156 


1.8 


64 


1.5e- 


05 


0.012 


3.1 


6.3e- 


-03 


0.354 


2.0 


2.7e- 


05 


0.022 


2.9 


l.le- 


-02 


0.597 


1.9 


128 


1.9e- 


06 


0.002 


3.0 


1.6e- 


-03 


0.088 


2.0 


3.7e- 


06 


0.003 


2.9 


3.4e- 


-03 


0.191 


1.6 



Table 3. Errors and rates of convergence for the test problem with polyno- 
mial solution using serendipity elements on asympotically affine meshes. 
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Uh\\L 2 




l|V(« 


- U h )\\ L 2 




n 


err. 


% 


rate 


err. 


% 


rate 


2 


5.0e- 


02 


4.066 




6.2e-01 


51.214 




4 


6.2e- 


03 


0.510 


3.0 


1.5c-01 


12.109 


2.1 


8 


7.6e- 


04 


0.062 


3.0 


3.6e-02 


2.948 


2.0 


16 


9.4e- 


05 


0.008 


3.0 


9.0e-03 


0.735 


2.0 


32 


1.2e- 


05 


0.001 


3.0 


2.2e-03 


0.183 


2.0 


64 


1.5e- 


06 


0.000 


3.0 


5.6e-04 


0.046 


2.0 


128 


1.9e- 


07 


0.000 


3.0 


1.4e-04 


0.012 


2.0 
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